%%%%%%%%%%%%%%%%%%
%% SETUP THE MATLAB PATHS
% make sure that fieldtrip and spm are not in your matlab path

% SET THE BELOW LINE TO THE OSL DIRECTORY:
%tilde='/home/mwoolrich';
tilde='/Users/woolrich';

osldir=[tilde '/homedir/matlab/osl1.3.1'];
addpath(osldir);
osl_startup(osldir);

%%%%%%%%%%%%%%%%%%
%% INITIALISE GLOBAL SETTINGS FOR THIS ANALYSIS
% This specifies where our data is stored, what the data filenames are, and
% a few parameters for preprocessing such as how we will epoch the data.

masksdir=[osldir '/std_masks']; % dir containing standard space masks
datadir=[tilde '/homedir/vols_data/george_data/alpha.oat']; % this is the directory the SPM files will be stored in

cmd = ['mkdir ' datadir]; unix(cmd); % make dir to put the results in

clear spm_files structural_files

% Alternatively you can leave fif_files empty and set up a list of SPM MEEG
% object files:
is=[1:3];
spm_files=[];
for i=1:length(is),    
    epoched_spm_files{i}=[datadir '/fsubject1' num2str(is(i)) '_spm_meeg.mat'];
end;

% structural files:
for i=1:length(is),
    structural_files{i}='';
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SETUP BEAMFORMER AND FIRST-LEVEL GLM OAT
% In this section we will do a wholebrain beamformer, followed by a trial-wise
% GLM that will correspond to a comparison of the ERFs for the different
% conditions.

%% SETUP THE OAT:
oat=osl_load_oat(datadir);

tmp=load([datadir '/fileRef']);
oat.subject_level.subjects_to_do=1:3;
oat.subject_level.session_index_list=tmp.ref.sessions;

oat.first_level.tf_method='hilbert';
oat.first_level.tf_num_freqs=1;
oat.first_level.tf_freq_range=[13 30];
oat.first_level.time_range=[-0.25 0.5];
oat.first_level.baseline_timespan=[-0.25 0];
oat.first_level.sessions_to_do=[1:9];
oat.first_level.tf_hilbert_freq_res=8;

oat.to_do=[0 0 0 1];

oat = osl_check_oat(oat);

oat = osl_run_oat(oat);

if(0),
    %%%%%%%%%%%%%%%%%%%
    %% OUTPUT SUBJECT'S NIFTII FILES
    % Having run the GLM on our source space data, we would like to inspect the
    % results for our single subject. 
    % We can do this by saving the contrast of parameter estimates (COPEs) and 
    % t-statistics for each of our contrasts to NIFTI images.

    S2=[];
    S2.oat=oat;
    S2.stats_fname=oat.first_level.results_fnames{1};
    S2.first_level_contrasts=[1,3]; % list of first level contrasts to output

    [statsdir times]=osl_save_nii_stats(S2);
end;

%%%%%%%%%%%%%%%%%%%
%% OUTPUT GROUP'S NIFTII FILES

S2=[];
S2.oat=oat;
S2.stats_fname=oat.group_level.results_fnames;
S2.first_level_contrasts=[3]; % list of first level contrasts to output

[statsdir,times]=osl_save_nii_stats(S2);
























%%%%%%%%%%%%%%%%%%%%%%%%
%% DO CLUSTER STATS
oat.group_level.cluster_stats_thresh=3;
oat.group_level.cluster_stats_nperms=300;
oat.to_do=[0 1 0];
oat = osl_run_oat(oat);

%% combine null distributions
    gstats_fname=oat.group_level.results_fnames;
    gstats=omt_load_stats(gstats_fname);

    distfile_save=[gstats_fname '_randomise_dir/permutations/dist' num2str(1) '.mat'];
    clusterstats=load(distfile_save);
    for i=2:oat.group_level.cluster_stats_nperms/100,
        distfile_save=[gstats_fname '_randomise_dir/permutations/dist' num2str(i) '.mat'];
        clusterstats_tmp=load(distfile_save);
        clusterstats.dist=[clusterstats.dist; clusterstats_tmp.dist];
    end;

    gstats.clusterstats{1,1}=clusterstats;
    omt_save_stats( gstats, gstats_fname );

    %% output nii files
    S2=[];
    S2.oat=oat;
    S2.stats_fname=oat.group_level.results_fnames;
    S2.first_level_contrasts=[1,3]; % list of first level contrasts to output
    S2.cluster_stats_fpr=0.05;
    statsdir=osl_save_nii_stats(S2);
   
   
    
    
%%%%%%%%%%%%%%%%%%%
%% OAT using ROI the whole way through

oat = osl_load_oat(oatdir); 

mask_fname=[OSLDIR '/std_masks/Right_Temporal_Occipital_Fusiform_Cortex'];
oat.group_level.time_range=[-0.1 0.3];
oat.group_level.time_average=0;
oat.group_level.time_smooth_std=0; % secs
oat.group_level.spatial_smooth_fwhm=0; % mm
oat.source_recon.mask_fname=[OSLDIR '/std_masks/Right_Temporal_Occipital_Fusiform_Cortex'];
oat.first_level.space_average=1;
oat.group_level.space_average=0;

oat.to_do=[1 1 1];

% run OAT
oat.group_level.use_robust_glm=0;
oat = osl_run_oat(oat);

gstats=osl_load_oat_results(oat,oat.group_level.results_fnames);

%%%%%%%%%%%%%%%%%%%
%% OAT using sensor space
S2=[];
S2.spm_files=spm_files;
S2.triggers={{'Motorbike'},{'Neutral face'},{'Happy face'},{'Fearful face'}};
S2.freq_range=[2 35]; % frequency range in Hz
S2.time_range=[-0.2 0.3];

% Xsummary is a parsimonious description of the design matrix.
% It contains values Xsummary{reg,cond}, where reg is a regressor no. and cond
% is a condition no. This will be used (by expanding the conditions over
% trials) to croat_settingse the (num_regressors x num_trials) design matrix:
Xsummary={};
Xsummary{1}=[1 0 0 0];Xsummary{2}=[0 1 0 0];Xsummary{3}=[0 0 1 0];Xsummary{4}=[0 0 0 1];
S2.design_matrix_summary=Xsummary;

% contrasts to be calculated:
S2.contrast={};
S2.contrast{1}=[3 0 0 0]'; % motorbikes
S2.contrast{2}=[0 1 1 1]'; % faces
S2.contrast{3}=[-3 1 1 1]'; % faces-motorbikes
S2.recon_method='none';
S2.modality_method='magplanar';
S2.time_range_stats=[-0.2 0.3];
S2.oat_name='/Users/woolrich/homedir/matlab/faces_group_data/faces_group_sensor';

oat = osl_setup_oat(S2);

%% run OAT
oat.to_do=[0 0 1];
oat.group_level.use_tstat=1;

oat = osl_run_oat(oat);

%% visualise using Fieldtrip
% note that this produces an interactive figure, with which you can:
% - draw around a set of sensors
% - click in the drawn box to produce a plot of the time series
% - on the time series plot you can draw a time window
% - and click in the window to create a topoplot averaged over that time
% window (which is itself interactive....!)

S2=[];
S2.oat=oat;
S2.stats_fname=oat.group_level.results_fnames;
S2.modality='mag';
S2.first_level_contrast=3;
S2.group_level_contrast=1;

% calculate t-stat using contrast of absolute value of parameter estimates
osl_stats_multiplotER(S2);


% SEE uploading_osl.txt in matlab dir


